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Abstract 

Using molecular dynamics we investigate the thermodynamics, dynamics and structure of 250 
diatomic molecules interacting by a core-softened potential. This system exhibits thermodynamics, 
dynamics and structural anomalies: a maximum in density-temperature plane at constante pressure 
and maximum and minimum points in the diffusivity and translational order parameter against 
density at constant temperature. Starting with very dense systems and decreasing density the 
mobility at low temperatures first increases, reach a maximum, then decreases, reach a minimum 
and finally increases. In the pressure-temperature phase diagram the line of maximum translational 
order parameter is located outside the line of diffusivity extrema that is enclosing the temperature 
of maximum density line. We compare our results with the monomeric system showing that the 
anisotropy due to the dumbbell leads to a much larger solid phase and to the appearance of a liquid 
crystal phase. 
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I. INTRODUCTION 



Several anomalies present in liquid water are believed to be due the formation and disrup- 
tion of hydrogen bonds. Specific angles and distances between water molecules are necessary 
to form such bonds, what leads to a dramatic competition between open low-density and 
closed high-density structures, depending on the thermodynamic state of the liquid. Al- 
though in water this competition actually acts in the second neighbors distances |lj], it could 
be mimetized as a competition between two preferred interparticle distances in a spherically 
symmetrical intermolecular interaction potential. Based on this mechanism, isotropic sys- 



tems where particle interact through core-softened potentials were considered 
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One of the biggest challenges for adopting these models to describe water and systems 
that exhibit the anomalies present in water is to reproduce some aspects of the anomalous 
behavior 
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291 ] . such as the existence of a temperature of maximum density (TMD), 



liquid-liquid phase transition or non-monotonic diffusion behavior with respect to isothermal 
pressure variation 
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The recently proposed core-softened shoulder potential [17], ]18|, ]19|, ]33j reproduces the 
density and diffusion anomalies. This potential can represent in an effective and orientation- 
averaged way the interaction between water pentamers characterized by the presence of two 
structures - one open and one closed - discussed above. Similarly, the thermodynamic and 
dynamic anomalies result from the competition between the two length scales associated 
with the open and closed structures. The open structure is favored by low pressures and 
the closed structure is favored by high pressures, but only becomes accessible at sufficiently 
high temperature. 

Simple pair potentials are particularly interesting because they are computatio nally 



cheaper than molecular models as well as amenable under analytica 
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381 ] is highly promis- 



For example, the recent perturbation theory developed by Zhou 
ing on attacking these sort of potentials. 

In s pite of core-softened potentials have been mainly used for modeling water 
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40, 4l|. 
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421 ]. many other materials present the so called water-like anomaly 
behaviour. For example, density anomaly was found experimentally in Se x Tei_ x ,[43] and 
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44j Liquid sulphur displays a sharp minimum in the density 45j , related to a 
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polymerization transition [46 
silica 



Waterlike anomalies were also found in simulations for 



50| silicon [5lJ and BeF 2 ,|47j. 
In this sense, it is reasonable to use core-softened potentials as the building blocks of 
a broader class of materials which we can classify as anomalous fluids. In principle this 
would allow us to fabricate anomalous liquids by imposing interparticle core-softened poten- 
tials. For instance, in principle one could build a polymer that diffuses faster under higher 
pressures what would be a quite interesting property for manufacturing plastic materials. 
But polymers are anisotropic systems and up to now the literature on core-softened poten- 
tials leading to water-like anomalies is restricted to spherical symmetric systems. Would 
the addition of an anisotropy in one direction to a system of particles interacting through 
a core-softened potential still maintain the anomalies present in the spherical symmetric 
system? Which new features the anisotropic system would exhibit? 

Here we address these questions by studying the pressure-temperature phase diagram of a 
model system made by rigid dimers. Each particle in the dimer interact with the particles in 
the other dimers through a core-softened potential. The results obtained using the dumbbell 
system are compared with pressure-temperature phase-diagram of the monomeric system 
interacting by the same core-softened potential employed in the dimeric case. 

Obviously anisotropic systems are not only computationally more complicated but also 
the addition of an extra degree of freedom yields a richer phase-diagram. For instance, 
diatomic particles interacting through a Lennard- Jones potential 52j exhibit a solid phase 
that occupies higher pressures and temperatures in the pressure-temperature phase diagram. 
In the case of the solid phase, the diatomic particles exhibit two close-packed arrangements 
instead of one observed in the monoatomic Lennard- Jones (53j. Therefore, we also expect the 
dimeric system interacting through two length scales potential might have a phase-diagram 
with a larger solid phase region in the pressure-temperature phase diagram than the one 



occupied by the solid phase in the monomeric 



system. Since in many cases the anomalies are 



close to the solid-liquid phase transition [2l|, |23[ the dumbbell model studied in this paper 



might have the anomalous region in the pressure-temperature phase diagram uncovered by 
the solid phases. We shall check if this is the case in this paper. 

The remaining of this manuscript goes as follows. In sec. [Ill the model is introduced and 
the simulation details are presented. In sec. II III the results are shown and conclusions are 
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found in sec. [TV] 



II. THE MODEL 



The model consists of N/2 dimeric molecules (dumbbells) formed each by two spherical 
particles of diameter a, linked rigidly in pairs with the distance A between their centers 
of mass as depicted in figure [TJ Each particle within a dimer interacts with all particles 
belonging to other dimers with the intermolecular continuous shoulder potential 18] given 
by 



U*(r) = 4 
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y of two length scales intermolecular interac 



r / V r , 
This potential can represent a whole fami 
tions, from a deep double wells potential 
the choice of the values of a, ro and c. 

This potential for the monomeric system was studied with a very small attractive region, 
with a — 5, To /cr = 0.7 and c = 1 so the liquid-gas unstable and metastable region would be 



331 ] to a repulsive shoulder [l3j, depending on 



avoided 
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19j. This potential has two length scales within a repulsive ramp followed by a 



very small attractive well. 

Here we explore the same interaction potential analyzed for the monomeric system, but 
for dumbbells. In particular we will assume X/a = 0.2. In order to study the equilibrium 
pressure-temperature phase diagram, we use molecular dynamics simulations to obtain the 
pressure as a function of temperature along isochores, diffusion constant as a function of 
density and temperature and the behavior of the structure as a function of temperature and 
pressure. 

We performed molecular dynamics simulations in the canonical ensemble using iV = 500 
particles (250 dimers) in a cubic box of volume V with periodic boundary conditions in 
the three directions, interacting with the intermolecular potential described above. The 
number density of the system is then p = N/V. The cutoff radius was set to 5.5er. Pressure, 
temperature, density, and diffusion are calculated in dimensionless units: 

T* = b*L 




r =r/a 



FIG. 1: Effective potential versus distance in reduced units. 
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In some state points we also carried out simulations with the same model but with 
1000 (500 dimers) and 2000 (1000 dimers) particles using the Large-scale Atomic/Molecular 
Massively Parallel Simulator (LAMMPS) J^] with essentially the same results. Further 
simulation details are discussed elsewhere [18]. 

Preliminary simulations showed that depending on the chosen temperature and density 
the system was in a fluid phase but became metastable with respect to the solid phase. 
In order to locate the phase boundary between the solid and the fluid phases two sets of 
simulations were carried out, one beginning with the molecules in a ordered crystal structure 
and the other beginning with the molecules in a random, liquid, starting structure obtained 
from previous equilibrium simulations. Thermodynamic and dynamic properties were cal- 
culated over 700 000 steps for the first set (and 900 000 steps for the second set), previously 
equilibrated over 200 000 (or 300 000) steps. The time step was 0.001 in reduced units and 



the time constant of the Berendsen thermostat [55j was 0.1 in reduced units. The internal 
bonds between the particles in each dimer remain fixed using the SHAKE jsfj algorithm, 
with a tolerance of 10~ 12 . 

The stability of the system was checked by analyzing the dependence of pressure on den- 
sity and also by visual analysis of the final structure, searching for cavitation. The structure 
of the system was characterized using the intermolecular radial distribution function, g(r) 
(RDF), which does not take into account the correlation between atoms belonging to the 
same molecule. The diffusion coefficient was calculated using the slope of the least square 
fit to the linear part of the mean square displacement, < r 2 {t) > (MSD), averaged over 
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different time origins. 

For analyzing the structure we define the structural anomaly region as the region where 
the translational order parameter t, given by 

t= IK, (3) 

Jo 

decreases upon increasing density. Here £ = rp 1 / 3 is the distance r in units of the mean 
interparticle separation p -1 / 3 , £ c is the cutoff distance set to half of the simulation box times 
p -1 / 3 , as in Ref. [lj]], <?(£) is the radial distribution function as a function of the (reduced) 
distance £ from a reference particle. For an ideal gas g = 1 and t = 0. In the crystal phase 
g ^ 1 over long distances and t is large. 

III. RESULTS 

A. The phase-diagram 

Fig. [2] shows (a) the pressure-temperature phase-diagram, (b) the radial distribution 
functions, (c) the mean square displacements and, finally, (d) snapshots of the system at some 
relevant thermodynamic state points. The pressure-temperature phase diagram, illustrated 
in Fig. E](a), displays at low temperatures a low density solid phase, a high density solid 
phase, a low density fluid phase, and a high density fluid phase. Near the boundaries of the 
high density solid phase, a liquid crystal-like (LCL) phase was also identified, as discussed 
below. 

At intermediate temperatures as the pressure is isothermally increased [following the 
arrow in the figure [2(a)] the system goes as follows: for low pressures the system is in the 
fluid phase, then as pressure increases it becomes solid, for higher pressures it becomes fluid 
again and then liquid crystal for even higher pressures. For very high pressures the system 
becomes solid and for even more higher pressures it becomes fluid again (not shown in the 
figure). 

The nature of the phases was investigated from three different ways: the radial distri- 
bution function, Fig. [2(b), the mean square displacement, Fig. [2(c), and the structural 
snapshots, Fig. [2(d). While for P* = 1.17 and 7.24 Fig. [2(d) shows ordered structures, 
for P* = 0.194 and 4.19 the structure is disordered, what could be typical for a liquid or a 
glass, depending on the mobility of particles in the system. According to the mean square 
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FIG. 2: (a) Reduced pressure versus reduced temperature phase-diagram showing the liquid and 
solid phases. The arrow at T* = 0.4 crosses the regions with P* = 0.194, 1.17, 4.19, and 7.24 
illustrated in the radial distribution function, (b) Radial distribution function versus reduced 
distance for T* = 0.4 and P* = 0.194, 1.17, 4.19, and 7.24. (c) Mean square displacement versus 
time for the four regions illustrated in (a): fluid, solid, fluid, and liquid crystal. With the slope of 
the MSD versus time for P* =0.194, 1.17, 4.19, 7.24 we calculate the diffusion coefficient that is 
zero for P* = 1.17. (d) Snapshots of the configurations for T* = 0.4 and the same pressures as in 
panel (a). 

displacement illustrated in Fig. [2](c), for P* = 0.194 and P* = 4.19 particles show high mo- 
bility, whereas for P* = 1.17 particles do not move. This would be already expected from 
the corresponding radial distribution function shown in Fig. [2(b). However, for P* = 7.24 
particles move almost as fast as in a liquid phase, but the corresponding radial distribution 
function and structural snapshots are rather solid-like. 

Which type of phase could be present in P* = 7.24? In order to answer this question a 
movie with the evolution of the configurations was made. This movie (available as supple- 
mentary material) shows that in this case particles actually move, but only in a string-like 
fashion, what characterizes a liquid crystal-like (LCL) phase. Besides the presence of this 
liquid crystal phase, the dumbbell system also exhibit a fluid phase at very high pressures. 
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B. The anomalies 

Similarly to the monomeric system previously studied [18j anomalies were also found in its 
dumbbells version, considered in this work. We focused in the three anomalies already found 
in the monomers case, i.e., the density, diffusion, and structural anomalies, as described 
below. 

(i) The density anomaly is the unusual expansion of the system upon cooling at constant 
pressure, a well-known effect which happens in water as discussed in the Introduction. In 
NPT-constant ensemble this anomaly is characterized by a maximum of the density along 
isobars in the density-temperature plane. Through thermodynamic relations 57j we are 
able to equivalently detect this anomaly by searching for a minimum of the pressure along 
isochores in the pressure-temperature phase-diagram. This was the technique used in this 
work since it is more suitable for the NVT ensemble. 

The density anomaly region for both (a) dimeric and (b) monomeric systems are shown in 
Fig. [3] as solid bold lines. Both dimeric and monomeric systems have the usual nose-shaped 



TMD line format, also found in molecular models for water 
isotropic potentials [22|, |39j, |40|, 141 ]. 
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3jj as well as for other 



For the dumbbells model the density anomaly region is located at higher pressures and 
temperatures and occupies a much larger region in the pressure-temperature phase-diagram 
than that observed for the the monomeric system. This can be evidenced by values for the 
maximum and minimum pressures and temperatures which enclose the TMD line for both 
systems, as shown in Table [H 

(ii) We also studied the molecular mobility calculating the diffusion coefficient using the 
the mean-square displacement averaged over different initial times, 



(Ar(t) 2 ) = ([r(t + t)-r(t )] 2 ). 



(4) 



Then the diffusion coefficient is obtained from the relation 



D = lim(Ar(t) 2 )/6t . 

t—>oo 



(5) 



Diffusion anomaly was also found in the dumbbells system. The usual procedure for 
detecting the dynamic anomaly region is to plot the diffusion coefficient versus density for 



fixed temperatures 
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571 ] as depicted in Fig. HI From this figure we see that for 
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FIG. 3: Pressure versus temperature phase diagram for (a) dimeric and (b) monomeric particles 
interacting through the potential illustrated in Fig. [J, and (c) the rescaled dimeric results, in which 
pressure and temperature from (a) are divided by 2 (see the text for more details). The results 
shown in panel (b) were adapted from Refs. [3] and The reentrant line with crosses in (a) 
represent the boundary between the fluid and the solid phases. The bold solid, dashed and dotted 
lines represent the TMD, D extrema, and the region of the structural anomaly respectively, in both 
(a) and (b) panels. 

a certain range of densities the diffusion coefficient anomalously increases under increasing 
density, changing its slope from negative to positive into this range. This determines a 
local maxima and minima for the D(p) plot. These extrema points can be mapped into 
a pressure-temperature plane, determining the lines of extrema in the diffusion coefficient 
inside which particles move faster under compression - or, equivalently, under increasing 
density. 

Fig. [3](a) illustrates the local diffusivity maxima and minima as dashed lines. Comparison 
between Fig. El^a) and Fig. E](b) indicates that the diffusion anomaly region for dumbbells 
occupies a larger region in the p-T plane than the diffusion anomaly region for the monomeric 
system. This can be evidenced by values for the maximum and minimum pressures and 
temperatures which enclose the diffusion anomalous region in Table [H 

(iii) We also study the behavior of the translational order parameter. We see from figure 
[5] that t decreases for increasing density in a certain range of densities, contrary of what 
expected for a normal fluid. The local maximum/minimum in the t(p) plot can be mapped 
into a pressure-temperature plane (as discussed above) giving the bounding lines for the 
region of structural anomaly in the P-T phase diagram, i.e., dotted lines in Fig. [3j 

Comparison between Fig. [3](a) and Fig. Mjo) shows that the structural anomaly region 
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FIG. 4: Diffusion coefficient versus density at fixed temperatures, which are T* =0.65, 0.75, 0.90, 
and 1.20 from bottom to top. 
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FIG. 5: Translational order parameter for the dumbbells system against density for fixed temper- 
atures, which are T* = 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, and 3.0 from top to bottom. 

is much broader for the dumbbells than for the monomers. It shall be stressed that the 
dumbbells system can be found in the liquid phase at pressures at high as P* = 17.6, which 
bounds the upper limit of the structural anomaly region. On the other hand, the lower limit 
for the structural anomaly region of the dumbbells model achieves pressures as low as 0.14 
in reduced units, even lower than for the monomers case (which is 0.32 in reduced units). 

In principle the differences between the monomeric and the dimeric system could be 
interpreted by assuming that a dimeric molecule is a monomeric with twice the potential e. 
In order to test if this would be the case, we rescaled pressure and temperature, T** = T*/2 
and P** = P*/2, and replot the diagram illustrated in Fig. [3] (a) obtaining Fig. [3] (c). 
Even in this case the dumbbell system shows remarkable differences when compared with 
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TABLE I: Maximum and minimum pressures and temperatures which bound (A) the TMD line 
(solid in Fig. [3]) in both systems, (B) dynamic anomaly region (bounded by dashed lines in Fig. 
[3]) for both systems and (C) the structural anomaly region (bounded by dotted lines in Fig. [3]) for 
both systems. 



ABC 
Dimers Monomers Dimers Monomers Dimers Monomers 



p* 


max. 


3.66 


0.90 


5.42 


0.94 


17.6 


0.94 


p* 


min. 


1.73 


0.55 


0.90 


0.43 


0.14 


0.43 


p* 


max. 


0.90 


0.26 


2.50 


0.40 


> 3.0 


1.00 


pjf 


min. 


0.65 


0.18 


0.70 


0.18 


< 0.70 


0.25 



the monomeric case. 

The addition of anisotropy changes the phase diagram of the fluid phase being the differ- 
ences significantly larger regarding the regions of the diffusion and structural anomalies. A 
proper choice of parameters, maybe a fine tuning of the distance A between the two particles 
of the same dimer, could lead to a specific range in the pressure-temperature phase space of 
the anomalies. 



IV. CONCLUSION 

In resume, we have addressed the question if just including anisotropy in a core-softened 
potential would lead to a modification in its anomalies and phases. For that purpose we 
have investigated a dimeric version of a previously studied monomeric model whose particles 



19] . 



were subject to a core-softened potential interaction [18|, 

We have found a much richer phase-diagram for the dimers than that one obtained with 
the monomers. While for the monomeric system the solid phase is at temperatures much 
lower than the temperatures of the TMD region, for the dumbbell case the solid phase 
occupies a much wider region in temperatures and terminates at the edge of the TMD 
region. This indicates that the anisotropy favors certain ordering of the system at low 
temperatures and low or intermediate pressures. In the case of the dimeric particles two 
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regions of solid phases separated by a liquid crystal region were observed, suggesting that 
the two length scales together with the anisotropic arrangement result in two solid phase 
densities. Even more surprising a very high pressure a liquid phase is also observed. 

The thermodynamic, dynamic, and structural anomalous regions maintain the same hi- 
erarchy in the dimeric case as observed in the monomeric system. The structural anomalous 
region occupies a wider region in the pressure temperature phase-diagram when compared 
with the diffusion anomalous region and this region occupies a wider region when compared 
with the density anomalous region. This hierarchy is also observed in water 
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. The 



dumbbell system, however, has pressures and temperatures in these anomalous regions much 
larger than the ones observed in the monomeric case. This behavior could be in principle 
explained mapping the dumbbell in a monomer with twice the potential interaction keeping 
the same interparticle distance. Since pressure and temperature scale with the interparticle 
potential, increasing the potential would lead to an increase in pressure and temperature. 

In resume, the addition of an anisotropy in a core-softened potential lead to a richer phase- 
diagram, but does not eliminates the thermodynamic, dynamic, and structural anomalies 
present in the isotropic system, but enlarges the solid state region. This result indicates 
that the use of dimers, trimers of even polymers of particles interacting with this isotropic 
potentials is a promissing way to design complex molecules which might lead to systems 
with the same anomalies present in water or even other kinds of anomalous behavior. 
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